Univariate Arima Models

Predicting COVID-19 Deaths (STAT 390)

1. Load Libraries, Seed, & Data

Code

# Load necessary libraries
library(forecast)
library(tseries)
library(tidyverse)
library(tidymodels)
library(kableExtra)
library(modeltime)
library(timetk)
library(cowplot)

# setting the seed
set.seed(123)

# Loading Data  ----------------------------------------------------------------
daily_deaths <- read_csv("data/new_daily_deaths.csv") %>% 
  select(daily_deaths)

east <- read_csv("data/east_daily.csv")

midwest <- read_csv("data/midwest_daily.csv")

south <- read_csv("data/south_daily.csv")

west <- read_csv("data/west_daily.csv")

2. Data Splitting

Code

# Data Splitting ---------------------------------------------------------------
east_splits <- time_series_split(east, initial = "1022 days", assess = "114 days")
east_train <- training(east_splits)
east_test <- testing(east_splits)

midwest_splits <- time_series_split(midwest, initial = "1022 days", assess = "114 days")
midwest_train <- training(midwest_splits)
midwest_test <- testing(midwest_splits)

south_splits <- time_series_split(south, initial = "1022 days", assess = "114 days")
south_train <- training(south_splits)
south_test <- testing(south_splits)

west_splits <- time_series_split(west, initial = "1022 days", assess = "114 days")
west_train <- training(west_splits)
west_test <- testing(west_splits)

## Convert Into Time Series ----------------------------------------------------
east_train_ts <- ts(east_train %>%
                      select(daily_deaths))
east_test_ts <- ts(east_test %>%
                     select(daily_deaths))
east_ts <- ts(east %>% 
                select(daily_deaths))

midwest_train_ts <- ts(midwest_train %>%
                      select(daily_deaths))
midwest_test_ts <- ts(midwest_test %>%
                     select(daily_deaths))
midwest_ts <- ts(midwest %>% 
                select(daily_deaths))

south_train_ts <- ts(south_train %>%
                    select(daily_deaths))
south_test_ts <- ts(south_test %>%
                   select(daily_deaths))
south_ts <- ts(south %>%
                 select(daily_deaths))

west_train_ts <- ts(west_train %>%
                   select(daily_deaths))
west_test_ts <- ts(west_test %>%
                  select(daily_deaths))
west_ts <- ts(west %>%
                select(daily_deaths))

## Plotting Training and Testing -----------------------------------------------
par(mfrow=c(2,2))
ggplot() +
  geom_line(data = east_test, aes(x = date, y = daily_deaths, color = "turquoise")) +
  geom_point(data = east_test, aes(x = date, y = daily_deaths, color = "turquoise"), size = .5) +
  geom_line(data = east_train, aes(x = date, y = daily_deaths, color = "orange")) +
  geom_point(data = east_train, aes(x = date, y = daily_deaths, color = "orange"), size = .5) +
  scale_color_manual(name="Data Set", values = c("turquoise", "orange"), labels = c("Training", "Testing")) +
  labs(y = "Daily Deaths",
       x = "Date",
       title = "East Region Training vs. Testing Data")

Code

ggplot() +
  geom_line(data = midwest_test, aes(x = date, y = daily_deaths, color = "turquoise")) +
  geom_point(data = midwest_test, aes(x = date, y = daily_deaths, color = "turquoise"), size = .5) +
  geom_line(data = midwest_train, aes(x = date, y = daily_deaths, color = "orange")) +
  geom_point(data = midwest_train, aes(x = date, y = daily_deaths, color = "orange"), size = .5) +
  scale_color_manual(name="Data Set", values = c("turquoise", "orange"), labels = c("Training", "Testing")) +
  labs(y = "Daily Deaths",
       x = "Date",
       title = "Midwest Region Training vs. Testing Data")

Code

ggplot() +
  geom_line(data = south_test, aes(x = date, y = daily_deaths, color = "turquoise")) +
  geom_point(data = south_test, aes(x = date, y = daily_deaths, color = "turquoise"), size = .5) +
  geom_line(data = south_train, aes(x = date, y = daily_deaths, color = "orange")) +
  geom_point(data = south_train, aes(x = date, y = daily_deaths, color = "orange"), size = .5) +
  scale_color_manual(name="Data Set", values = c("turquoise", "orange"), labels = c("Training", "Testing")) +
  labs(y = "Daily Deaths",
       x = "Date",
       title = "South Region Training vs. Testing Data")

Code

ggplot() +
  geom_line(data = west_test, aes(x = date, y = daily_deaths, color = "turquoise")) +
  geom_point(data = west_test, aes(x = date, y = daily_deaths, color = "turquoise"), size = .5) +
  geom_line(data = west_train, aes(x = date, y = daily_deaths, color = "orange")) +
  geom_point(data = west_train, aes(x = date, y = daily_deaths, color = "orange"), size = .5) +
  scale_color_manual(name="Data Set", values = c("turquoise", "orange"), labels = c("Training", "Testing")) +
  labs(y = "Daily Deaths",
       x = "Date",
       title = "West Region Training vs. Testing Data")

3. Testing for Stationarity

Code

### Using Augmented Dickey-Fuller Test -----------------------------------------
adf_east <- adf.test(east_ts)
print(adf_east) # p-value is 0.01 < 0.05, implying it is stationary

    Augmented Dickey-Fuller Test

data:  east_ts
Dickey-Fuller = -4.5566, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary

Code

adf_midwest <- adf.test(midwest_ts)
print(adf_midwest) # p-value is 0.4919 > 0.05, implying it is NOT stationary

    Augmented Dickey-Fuller Test

data:  midwest_ts
Dickey-Fuller = -2.5984, Lag order = 10, p-value = 0.325
alternative hypothesis: stationary

Code

midwest_ts <- diff(midwest_ts)
adf_midwest <- adf.test(midwest_ts)
print(adf_midwest) # p-value is 0.01 < 0.05, implying it is stationary with 1 diff

    Augmented Dickey-Fuller Test

data:  midwest_ts
Dickey-Fuller = -13.43, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary

Code

adf_south <- adf.test(south_ts)
print(adf_south) # p-value is 0.3018 > 0.05, implying it is NOT stationary

    Augmented Dickey-Fuller Test

data:  south_ts
Dickey-Fuller = -2.6533, Lag order = 10, p-value = 0.3018
alternative hypothesis: stationary

Code

south_ts <- diff(south_ts)
adf_south <- adf.test(south_ts)
print(adf_south) # p-value is 0.01 < 0.05, implying it is stationary with 1 diff

    Augmented Dickey-Fuller Test

data:  south_ts
Dickey-Fuller = -11.36, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary

Code

adf_west <- adf.test(west_ts)
print(adf_west) # p-value is 0.4098 > 0.05, implying it is NOT stationary

    Augmented Dickey-Fuller Test

data:  west_ts
Dickey-Fuller = -2.3981, Lag order = 10, p-value = 0.4098
alternative hypothesis: stationary

Code

west_ts <- diff(west_ts)
adf_west <- adf.test(west_ts)
print(adf_west)  # p-value is 0.01 < 0.05, implying it is stationary with 1 diff

    Augmented Dickey-Fuller Test

data:  west_ts
Dickey-Fuller = -12.247, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary

Code

### Making Testing/Training Stationary -----------------------------------------
midwest_train_ts <- diff(midwest_train_ts)
midwest_test_ts <- diff(midwest_test_ts)

south_train_ts <- diff(south_train_ts)
south_test_ts <- diff(south_test_ts)

west_train_ts <- diff(west_train_ts)
west_test_ts <- diff(west_test_ts)

4. ACF & PACF Plots

Code

## ACF & PACF Plots ------------------------------------------------------------
par(mfrow=c(3, 1))
acf(east_ts, main = "East")
pacf(east_ts, main = "East")
plot(east_ts, main = "Daily Deaths: East")

Code

acf(midwest_ts, main = "Midwest")
pacf(midwest_ts, main = "Midwest")
plot(midwest_ts, main = "Daily Deaths: Midwest")

Code

acf(south_ts, main = "South")
pacf(south_ts, main = "South")
plot(south_ts, main = "Daily Deaths: South")

Code

acf(west_ts, main = "West")
pacf(west_ts, main = "West")
plot(south_ts, main = "Daily Deaths: South")

Code

plot(west_ts, main = "Daily Deaths: West")

5. Time Series Decomposition

Code

# total daily death decomposition
daily_ts <- ts(daily_deaths, frequency = 365)

decomposed_ts <- decompose(daily_ts)

# Plot the original time series and the decomposed components
par(mfrow=c(2,1))
plot(daily_ts, main = "Total Original Time Series", ylab = "Value", col = "blue")
plot(decomposed_ts$trend, main = "Trend Component", ylab = "Value", col = "red")

Code

plot(decomposed_ts$seasonal, main = "Seasonal Component", ylab = "Value", col = "green")
plot(decomposed_ts$random, main = "Residual Component", ylab = "Value", col = "purple")

Code

# East Decomposition -----------------------------------------------------------
east_daily <- east %>% 
  select(daily_deaths)

east_ts <- ts(east_daily, frequency = 365)

decomposed_east <- decompose(east_ts)

# Plot the original time series and the decomposed components
par(mfrow=c(2,1))
plot(east_ts, main = "East Original Time Series", ylab = "Value", col = "blue")
plot(decomposed_east$trend, main = "Trend Component", ylab = "Value", col = "red")

Code

plot(decomposed_east$seasonal, main = "Seasonal Component", ylab = "Value", col = "green")
plot(decomposed_east$random, main = "Residual Component", ylab = "Value", col = "purple")

Code

## Midwest Decomposition ----------------------------------------------------------
midwest_daily <- midwest %>% 
  select(daily_deaths)

midwest_ts <- ts(midwest_daily, frequency = 365)

decomposed_midwest <- decompose(midwest_ts)

# Plot the original time series and the decomposed components
par(mfrow=c(2,1))
plot(midwest_ts, main = "Midwest Original Time Series", ylab = "Value", col = "blue")
plot(decomposed_midwest$trend, main = "Trend Component", ylab = "Value", col = "red")

Code

plot(decomposed_midwest$seasonal, main = "Seasonal Component", ylab = "Value", col = "green")
plot(decomposed_midwest$random, main = "Residual Component", ylab = "Value", col = "purple")

Code

## South Decomposition ----------------------------------------------------------
south_daily <- south %>% 
  select(daily_deaths)

south_ts <- ts(south_daily, frequency = 365)

decomposed_south <- decompose(south_ts)

# Plot the original time series and the decomposed components
par(mfrow=c(2,1))
plot(south_ts, main = "South Original Time Series", ylab = "Value", col = "blue")
plot(decomposed_south$trend, main = "Trend Component", ylab = "Value", col = "red")

Code

plot(decomposed_south$seasonal, main = "Seasonal Component", ylab = "Value", col = "green")
plot(decomposed_south$random, main = "Residual Component", ylab = "Value", col = "purple")

Code

## West Decomposition ----------------------------------------------------------
west_daily <- west %>% 
  select(daily_deaths)

west_ts <- ts(west_daily, frequency = 365)

decomposed_west <- decompose(west_ts)

# Plot the original time series and the decomposed components
par(mfrow=c(2,1))
plot(west_ts, main = "West Original Time Series", ylab = "Value", col = "blue")
plot(decomposed_west$trend, main = "Trend Component", ylab = "Value", col = "red")

Code

plot(decomposed_west$seasonal, main = "Seasonal Component", ylab = "Value", col = "green")
plot(decomposed_west$random, main = "Residual Component", ylab = "Value", col = "purple")

6. ARIMA Model

6.1 Initial Model

Code

# Building ARIMA Model ---------------------------------------------------------
east_arima <- arima(east_train_ts, order=c(1,0,1))
midwest_arima <- arima(midwest_train_ts, order=c(1,1,1))
south_arima <- arima(south_train_ts, order=c(1,1,1))
west_arima <- arima(west_train_ts, order=c(1,1,1))

## Summary & Forecast ---------------------------------------------------------
summary(east_arima)

Call:
arima(x = east_train_ts, order = c(1, 0, 1))

Coefficients:
         ar1      ma1  intercept
      0.9176  -0.0388   197.3738
s.e.  0.0147   0.0493    39.7998

sigma^2 estimated as 12153:  log likelihood = -6257.17,  aic = 12522.33

Training set error measures:
                    ME     RMSE      MAE  MPE MAPE      MASE        ACF1
Training set 0.1151898 110.2407 64.66739 -Inf  Inf 0.9918341 0.008515868

Code

summary(midwest_arima)

Call:
arima(x = midwest_train_ts, order = c(1, 1, 1))

Coefficients:
          ar1      ma1
      -0.3477  -1.0000
s.e.   0.0293   0.0026

sigma^2 estimated as 50173:  log likelihood = -6970.99,  aic = 13947.98

Training set error measures:
                     ME     RMSE      MAE MPE MAPE      MASE        ACF1
Training set -0.8748943 223.8834 138.9352 NaN  Inf 0.6215031 -0.06581815

Code

summary(south_arima)

Call:
arima(x = south_train_ts, order = c(1, 1, 1))

Coefficients:
          ar1      ma1
      -0.1727  -1.0000
s.e.   0.0308   0.0026

sigma^2 estimated as 106737:  log likelihood = -7355.8,  aic = 14717.6

Training set error measures:
                    ME    RMSE      MAE MPE MAPE      MASE        ACF1
Training set -1.420592 326.547 203.7282 NaN  Inf 0.6344327 -0.04113314

Code

summary(west_arima)

Call:
arima(x = west_train_ts, order = c(1, 1, 1))

Coefficients:
          ar1      ma1
      -0.1698  -1.0000
s.e.   0.0309   0.0026

sigma^2 estimated as 23955:  log likelihood = -6593.77,  aic = 13193.54

Training set error measures:
                     ME     RMSE      MAE MPE MAPE      MASE        ACF1
Training set -0.6839532 154.6997 101.0808 NaN  Inf 0.6424269 -0.02682684

Code

east_values <- forecast(east_arima, h = length(east_test_ts))
midwest_values <- forecast(midwest_arima, h = length(midwest_test_ts))
south_values <- forecast(south_arima, h = length(south_test_ts))
west_values <- forecast(west_arima, h = length(west_test_ts))

6.2 Accuracy

Code

## Evaluating Accuracy ---------------------------------------------------------
east_initial_arima_acc <- forecast::accuracy(east_values)

east_initial_arima_acc %>% 
  kbl(caption = "<center>Initial East Accuracy<center>") %>% 
  kable_styling()

Initial East Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.1151898	110.2407	64.66739	-Inf	Inf	0.9918341	0.0085159

Code

### East Region
east_test_plot <- as.data.frame(east_test_ts) %>% 
  mutate(num = seq(1023, (1023+113)))

east_initial_arima_plot <- forecast::autoplot(east_values) +
  geom_line(data = east_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "East Region Initial ARIMA",
       y = "Daily Deaths")

east_initial_arima_plot

Code

### Midwest Region

midwest_initial_arima_acc <- forecast::accuracy(midwest_values)

midwest_initial_arima_acc %>% 
  kbl(caption = "<center>Initial Midwest Accuracy<center>") %>% 
  kable_styling()

Initial Midwest Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	-0.8748943	223.8834	138.9352	NaN	Inf	0.6215031	-0.0658182

Code

midwest_test_plot <- as.data.frame(midwest_test_ts) %>% 
  mutate(num = seq(1023, (1023+112)))

midwest_initial_arima_plot <- forecast::autoplot(midwest_values) +
  geom_line(data = midwest_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "Midwest Region Initial ARIMA",
       y = "Daily Deaths")

midwest_initial_arima_plot

Code

### South Region

south_initial_arima_acc <- forecast::accuracy(south_values)

south_initial_arima_acc %>% 
  kbl(caption = "<center>Initial South Accuracy<center>") %>% 
  kable_styling()

Initial South Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	-1.420592	326.547	203.7282	NaN	Inf	0.6344327	-0.0411331

Code

south_test_plot <- as.data.frame(south_test_ts) %>% 
  mutate(num = seq(1023, (1023+112)))

south_initial_arima_plot <- forecast::autoplot(south_values) +
  geom_line(data = south_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "South Region Initial ARIMA",
       y = "Daily Deaths")

south_initial_arima_plot 

Code

### West Region

west_initial_arima_acc <- forecast::accuracy(west_values)

west_initial_arima_acc %>% 
  kbl(caption = "<center>Initial West Accuracy<center>") %>% 
  kable_styling()

Initial West Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	-0.6839532	154.6997	101.0808	NaN	Inf	0.6424269	-0.0268268

Code

west_test_plot <- as.data.frame(west_test_ts) %>% 
  mutate(num = seq(1023, (1023+112)))

west_initial_arima_plot <- forecast::autoplot(west_values) +
  geom_line(data = west_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "West Region Initial ARIMA",
       y = "Daily Deaths") 

west_initial_arima_plot

6.3 Grid Search

Code

# Grid Search for p and q ------------------------------------------------------

p_values <- 0:5  # AR order
d_values <- 0:0  # I order
q_values <- 0:5  # MA order

## East - Perform grid search for ARIMA parameters  ----------------------------
best_east_model <- NULL
best_east_aic <- Inf

for (p in p_values) {
  for (d in d_values) {
    for (q in q_values) {
      # Fit ARIMA model
      arima_model <- tryCatch(
        {
          fit <- Arima(east_ts, order = c(p, d, q))
          fit
        },
        error = function(e) {
          NULL
        }
      )
      
      # Check if ARIMA model was successfully fitted
      if (!is.null(arima_model)) {
        # Check AIC criterion
        if (AIC(arima_model) < best_east_aic) {
          best_east_model <- arima_model
          best_east_aic <- AIC(arima_model)
          best_east_order <- c(p, d, q)
        }
      }
    }
  }
}

Code

## Midwest - Perform grid search for ARIMA parameters  -------------------------
p_values <- 0:5  # AR order
d_values <- 0:1  # I order
q_values <- 0:5  # MA order

best_midwest_model <- NULL
best_midwest_aic <- Inf

for (p in p_values) {
  for (d in d_values) {
    for (q in q_values) {
      # Fit ARIMA model
      arima_model <- tryCatch(
        {
          fit <- Arima(midwest_ts, order = c(p, d, q))
          fit
        },
        error = function(e) {
          NULL
        }
      )
      
      # Check if ARIMA model was successfully fitted
      if (!is.null(arima_model)) {
        # Check AIC criterion
        if (AIC(arima_model) < best_midwest_aic) {
          best_midwest_model <- arima_model
          best_midwest_aic <- AIC(arima_model)
          best_midwest_order <- c(p, d, q)
        }
      }
    }
  }
}

Code

## South - Perform grid search for ARIMA parameters  ---------------------------
p_values <- 0:5  # AR order
d_values <- 0:1  # I order
q_values <- 0:5  # MA order

best_south_model <- NULL
best_south_aic <- Inf

for (p in p_values) {
  for (d in d_values) {
    for (q in q_values) {
      # Fit ARIMA model
      arima_model <- tryCatch(
        {
          fit <- Arima(south_ts, order = c(p, d, q))
          fit
        },
        error = function(e) {
          NULL
        }
      )
      
      # Check if ARIMA model was successfully fitted
      if (!is.null(arima_model)) {
        # Check AIC criterion
        if (AIC(arima_model) < best_south_aic) {
          best_south_model <- arima_model
          best_south_aic <- AIC(arima_model)
          best_south_order <- c(p, d, q)
        }
      }
    }
  }
}

Code

## West - Perform grid search for ARIMA parameters  ----------------------------
p_values <- 0:5  # AR order
d_values <- 0:1  # I order
q_values <- 0:5  # MA order

best_west_model <- NULL
best_west_aic <- Inf

for (p in p_values) {
  for (d in d_values) {
    for (q in q_values) {
      # Fit ARIMA model
      arima_model <- tryCatch(
        {
          fit <- Arima(west_ts, order = c(p, d, q))
          fit
        },
        error = function(e) {
          NULL
        }
      )
      
      # Check if ARIMA model was successfully fitted
      if (!is.null(arima_model)) {
        # Check AIC criterion
        if (AIC(arima_model) < best_west_aic) {
          best_west_model <- arima_model
          best_west_aic <- AIC(arima_model)
          best_west_order <- c(p, d, q)
        }
      }
    }
  }
}

Code

# Saving results
save(best_east_order, best_east_model,
     best_south_order, best_south_model,
     best_midwest_order, best_midwest_model,
     best_west_order, best_west_model, 
     file = "results/arima_grid.rda")

Best East Model

Code

# load results
load("results/arima_grid.rda")
# Display the best ARIMA model and its parameters
print(paste(best_east_order))

[1] "5" "0" "5"

Code

print(best_east_model)

Series: east_ts 
ARIMA(5,0,5) with non-zero mean 

Coefficients:
         ar1      ar2     ar3     ar4      ar5      ma1     ma2     ma3
      2.5325  -2.4820  0.7383  0.6063  -0.4009  -2.1286  1.6389  0.0590
s.e.  0.1056   0.3215  0.4321  0.3068   0.0949   0.0975  0.2582  0.3121
          ma4     ma5      mean
      -0.8826  0.4627  184.1912
s.e.   0.1987  0.0577   60.5578

sigma^2 = 6735:  log likelihood = -6616.45
AIC=13256.9   AICc=13257.17   BIC=13317.32

Best Midwest Model

Code

print(paste(best_midwest_order))

[1] "5" "1" "5"

Code

print(best_midwest_model)

Series: midwest_ts 
ARIMA(5,1,5) 

Coefficients:
         ar1      ar2      ar3      ar4      ar5      ma1     ma2      ma3
      0.2188  -0.9672  -0.0252  -0.5260  -0.5481  -1.1980  1.4606  -1.2997
s.e.  0.0353   0.0286   0.0485   0.0288   0.0324   0.0399  0.0491   0.0522
         ma4      ma5
      1.0143  -0.4534
s.e.  0.0483   0.0337

sigma^2 = 15644:  log likelihood = -7089.78
AIC=14201.56   AICc=14201.8   BIC=14256.94

Best South Model

Code

print(paste(best_south_order))

[1] "5" "1" "5"

Code

print(best_south_model)

Series: south_ts 
ARIMA(5,1,5) 

Coefficients:
         ar1      ar2     ar3      ar4      ar5      ma1     ma2      ma3
      0.4802  -1.1648  0.3304  -0.7153  -0.3031  -1.3248  1.5824  -1.3808
s.e.  0.0544   0.0456  0.0758   0.0423   0.0488   0.0519  0.0749   0.0844
         ma4      ma5
      1.0487  -0.4103
s.e.  0.0652   0.0413

sigma^2 = 38413:  log likelihood = -7598.92
AIC=15219.84   AICc=15220.08   BIC=15275.22

Best West Model

Code

print(paste(best_west_order))

[1] "5" "1" "5"

Code

print(best_west_model)

Series: west_ts 
ARIMA(5,1,5) 

Coefficients:
         ar1      ar2     ar3      ar4      ar5      ma1     ma2      ma3
      0.4629  -1.1476  0.3025  -0.7075  -0.3255  -1.3275  1.6081  -1.4132
s.e.  0.0516   0.0452  0.0734   0.0403   0.0465   0.0492  0.0746   0.0849
         ma4      ma5
      1.1110  -0.4753
s.e.  0.0595   0.0368

sigma^2 = 8126:  log likelihood = -6717.71
AIC=13457.43   AICc=13457.66   BIC=13512.8

6.4 Rebuilding ARIMA Model

Code

# Rebuilding ARIMA Models ------------------------------------------------------
east_arima <- arima(east_train_ts, order=c(5,0,5))
midwest_arima <- arima(midwest_train_ts, order=c(5,1,5))
south_arima <- arima(south_train_ts, order=c(5,1,5))
west_arima <- arima(west_train_ts, order=c(5,1,5), method = "ML")


east_values <- forecast(east_arima, h = length(east_test_ts))
midwest_values <- forecast(midwest_arima, h = length(midwest_test_ts))
south_values <- forecast(south_arima, h = length(south_test_ts))
west_values <- forecast(west_arima, h = length(west_test_ts))

## Evaluating Accuracy & Plots -------------------------------------------------
east_tuned_arima_acc <- forecast::accuracy(east_values)

east_tuned_arima_acc %>% 
  kbl(caption = "<center>Tuned East Accuracy<center>") %>% 
  kable_styling()

Tuned East Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.1095068	80.35755	47.18823	NaN	Inf	0.723748	-0.0254917

Code

### East Region
east_tuned_arima_plot <- forecast::autoplot(east_values) +
  geom_line(data = east_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "East Region Tuned ARIMA",
       y = "Daily Deaths")

east_tuned_arima_plot

Code

### Midwest Region

midwest_tuned_arima_acc <- forecast::accuracy(midwest_values)

midwest_tuned_arima_acc %>% 
  kbl(caption = "<center>Tuned Midwest Accuracy<center>") %>% 
  kable_styling()

Tuned Midwest Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	-3.635509	152.8728	95.00769	NaN	Inf	0.4250009	-0.0128029

Code

midwest_tuned_arima_plot <- forecast::autoplot(midwest_values) +
  geom_line(data = midwest_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "Midwest Region Tuned ARIMA",
       y = "Daily Deaths")

midwest_tuned_arima_plot

Code

### South Region

south_tuned_arima_acc <- forecast::accuracy(south_values)

south_tuned_arima_acc %>% 
  kbl(caption = "<center>Tuned South Accuracy<center>") %>% 
  kable_styling()

Tuned South Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	-4.199063	200.7617	120.224	NaN	Inf	0.3743912	-0.0666043

Code

south_tuned_arima_plot <- forecast::autoplot(south_values) +
  geom_line(data = south_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "South Region Tuned ARIMA",
       y = "Daily Deaths")

south_tuned_arima_plot 

Code

### West Region
west_tuned_arima_acc <- forecast::accuracy(west_values)

west_tuned_arima_acc %>% 
  kbl(caption = "<center>Tuned West Accuracy<center>") %>% 
  kable_styling()

Tuned West Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	-0.0352798	99.27405	63.93076	NaN	Inf	0.4063168	-0.0093102

Code

west_tuned_arima_plot <- forecast::autoplot(west_values) +
  geom_line(data = west_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "West Region Tuned ARIMA",
       y = "Daily Deaths") 

west_tuned_arima_plot

Code

# saving results
save(east_tuned_arima_acc, east_initial_arima_acc, east_initial_arima_plot, east_tuned_arima_plot,
     midwest_tuned_arima_acc, midwest_initial_arima_acc, midwest_initial_arima_plot, midwest_tuned_arima_plot,
     south_tuned_arima_acc, south_initial_arima_acc, south_initial_arima_plot, south_tuned_arima_plot,
     west_tuned_arima_acc, west_initial_arima_acc, west_initial_arima_plot, west_tuned_arima_plot,
     file = "results/uni_arima.rda"
     )

7. SARIMA Model

7.1 Initial SARIMA Model

Code

# Building SARIMA model  ------------------------------------------------------
east_sarima <- arima(east_train_ts, seasonal = list(order = c(1,1,1)))
midwest_sarima <- arima(midwest_train_ts, seasonal = list(order = c(1,1,1)))
south_sarima <- arima(south_train_ts, seasonal = list(order = c(1,1,1)))
west_sarima <- arima(west_train_ts, seasonal = list(order = c(1,1,1)))

# Summary of the SARIMA model
summary(east_sarima)

Call:
arima(x = east_train_ts, seasonal = list(order = c(1, 1, 1)))

Coefficients:
        sar1     sma1
      0.4706  -0.8016
s.e.  0.0382   0.0211

sigma^2 estimated as 11156:  log likelihood = -6206.64,  aic = 12419.27

Training set error measures:
                    ME     RMSE     MAE  MPE MAPE      MASE      ACF1
Training set 0.1606608 105.5715 61.2767 -Inf  Inf 0.9398295 0.1140208

Code

summary(midwest_sarima)

Call:
arima(x = midwest_train_ts, seasonal = list(order = c(1, 1, 1)))

Coefficients:
         sar1     sma1
      -0.3477  -1.0000
s.e.   0.0293   0.0026

sigma^2 estimated as 50173:  log likelihood = -6970.99,  aic = 13947.98

Training set error measures:
                     ME     RMSE      MAE MPE MAPE      MASE        ACF1
Training set -0.8748943 223.8834 138.9352 NaN  Inf 0.6215031 -0.06581815

Code

summary(south_sarima)

Call:
arima(x = south_train_ts, seasonal = list(order = c(1, 1, 1)))

Coefficients:
         sar1     sma1
      -0.1727  -1.0000
s.e.   0.0308   0.0026

sigma^2 estimated as 106737:  log likelihood = -7355.8,  aic = 14717.6

Training set error measures:
                    ME    RMSE      MAE MPE MAPE      MASE        ACF1
Training set -1.420592 326.547 203.7282 NaN  Inf 0.6344327 -0.04113314

Code

summary(west_sarima)

Call:
arima(x = west_train_ts, seasonal = list(order = c(1, 1, 1)))

Coefficients:
         sar1     sma1
      -0.1698  -1.0000
s.e.   0.0309   0.0026

sigma^2 estimated as 23955:  log likelihood = -6593.77,  aic = 13193.54

Training set error measures:
                     ME     RMSE      MAE MPE MAPE      MASE        ACF1
Training set -0.6839532 154.6997 101.0808 NaN  Inf 0.6424269 -0.02682684

Code

## Forecasting with the SARIMA model -------------------------------------------
east_values <- forecast(east_sarima, h = length(east_test_ts))
midwest_values <- forecast(midwest_sarima, h = length(midwest_test_ts))
south_values <- forecast(south_sarima, h = length(south_test_ts))
west_values <- forecast(west_sarima, h = length(west_test_ts))

7.2 Accuracy

Code

## Plotting & Accuracy

## East Region
east_initial_sarima_acc <- forecast::accuracy(east_values)

east_initial_sarima_acc %>% 
  kbl(caption = "<center>Initial SARIMA East Accuracy<center>") %>% 
  kable_styling()

Initial SARIMA East Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.1606608	105.5715	61.2767	-Inf	Inf	0.9398295	0.1140208

Code

east_initial_sarima_plot <- forecast::autoplot(east_values) +
  geom_line(data = east_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "East Region Initial SARIMA",
       y = "Daily Deaths")

east_initial_sarima_plot

Code

## Midwest Region
midwest_initial_sarima_acc <- forecast::accuracy(midwest_values)

midwest_initial_sarima_acc %>% 
  kbl(caption = "<center>Initial SARIMA Midwest Accuracy<center>") %>% 
  kable_styling()

Initial SARIMA Midwest Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	-0.8748943	223.8834	138.9352	NaN	Inf	0.6215031	-0.0658182

Code

midwest_initial_sarima_plot <- forecast::autoplot(midwest_values) +
  geom_line(data = midwest_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "Midwest Region Initial SARIMA",
       y = "Daily Deaths")

midwest_initial_sarima_plot

Code

## South Region
south_initial_sarima_acc <- forecast::accuracy(south_values)

south_initial_sarima_acc %>% 
  kbl(caption = "<center>Initial SARIMA South Accuracy<center>") %>% 
  kable_styling()

Initial SARIMA South Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	-1.420592	326.547	203.7282	NaN	Inf	0.6344327	-0.0411331

Code

south_initial_sarima_plot <- forecast::autoplot(south_values) +
  geom_line(data = south_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "South Region Initial SARIMA",
       y = "Daily Deaths")

south_initial_sarima_plot

Code

## West Region
west_initial_sarima_acc <- forecast::accuracy(west_values)

west_initial_sarima_acc %>% 
  kbl(caption = "<center>Initial SARIMA West Accuracy<center>") %>% 
  kable_styling()

Initial SARIMA West Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	-0.6839532	154.6997	101.0808	NaN	Inf	0.6424269	-0.0268268

Code

west_initial_sarima_plot <- forecast::autoplot(west_values) +
  geom_line(data = west_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "West Region Initial SARIMA",
       y = "Daily Deaths")

west_initial_sarima_plot

7.3 Grid Search

Code

## Define the parameter grids --------------------------------------------------
p_grid <- 0:5  # AR parameter
d_grid <- 0:0  # Differencing parameter
q_grid <- 0:5  # MA parameter
P_grid <- 0:1  # Seasonal AR parameter
D_grid <- 0:1  # Seasonal differencing parameter
Q_grid <- 0:1  # Seasonal MA parameter
s <- 12        # Seasonal period



## East Grid Search ------------------------------------------------------------

# Create East Region
east_train <- east_train %>% 
  select(daily_deaths)

# Create an empty data frame to store the results
east_results <- data.frame(order = character(),
                      seasonal = character(),
                      AIC = numeric(),
                      stringsAsFactors = FALSE)

# Perform grid search
for (p in p_grid) {
  for (d in d_grid) {
    for (q in q_grid) {
      for (P in P_grid) {
        for (D in D_grid) {
          for (Q in Q_grid) {
            order <- c(p, d, q)
            seasonal <- c(P, D, Q, s)
            model <- tryCatch(
              {
                fit <- arima(east_train, order = order, seasonal = seasonal)
                AIC_val <- AIC(fit)
                # Append the results to the data frame
                east_results <- rbind(east_results, data.frame(order = paste(order, collapse = ", "),
                                                     seasonal = paste(seasonal, collapse = ", "),
                                                     AIC = AIC_val))
              },
              error = function(e) NULL
            )
          }
        }
      }
    }
  }
}

d_grid <- 0:1  # Differencing parameter

## Midwest Grid Search ------------------------------------------------------------

# Create Midwest Region
midwest_train <- midwest_train %>% 
  select(daily_deaths)

# Create an empty data frame to store the results
midwest_results <- data.frame(order = character(),
                           seasonal = character(),
                           AIC = numeric(),
                           stringsAsFactors = FALSE)

# Perform grid search
for (p in p_grid) {
  for (d in d_grid) {
    for (q in q_grid) {
      for (P in P_grid) {
        for (D in D_grid) {
          for (Q in Q_grid) {
            order <- c(p, d, q)
            seasonal <- c(P, D, Q, s)
            model <- tryCatch(
              {
                fit <- arima(midwest_train, order = order, seasonal = seasonal)
                AIC_val <- AIC(fit)
                # Append the results to the data frame
                midwest_results <- rbind(midwest_results, data.frame(order = paste(order, collapse = ", "),
                                                               seasonal = paste(seasonal, collapse = ", "),
                                                               AIC = AIC_val))
              },
              error = function(e) NULL
            )
          }
        }
      }
    }
  }
}

## South Grid Search ------------------------------------------------------------

# Create South Region
south_train <- south_train %>% 
  select(daily_deaths)

# Create an empty data frame to store the results
south_results <- data.frame(order = character(),
                           seasonal = character(),
                           AIC = numeric(),
                           stringsAsFactors = FALSE)

# Perform grid search
for (p in p_grid) {
  for (d in d_grid) {
    for (q in q_grid) {
      for (P in P_grid) {
        for (D in D_grid) {
          for (Q in Q_grid) {
            order <- c(p, d, q)
            seasonal <- c(P, D, Q, s)
            model <- tryCatch(
              {
                fit <- arima(south_train, order = order, seasonal = seasonal)
                AIC_val <- AIC(fit)
                # Append the results to the data frame
                south_results <- rbind(south_results, data.frame(order = paste(order, collapse = ", "),
                                                               seasonal = paste(seasonal, collapse = ", "),
                                                               AIC = AIC_val))
              },
              error = function(e) NULL
            )
          }
        }
      }
    }
  }
}


## West Grid Search ------------------------------------------------------------

# Create West Region
west_train <- west_train %>% 
  select(daily_deaths)

# Create an empty data frame to store the results
west_results <- data.frame(order = character(),
                           seasonal = character(),
                           AIC = numeric(),
                           stringsAsFactors = FALSE)

# Perform grid search
for (p in p_grid) {
  for (d in d_grid) {
    for (q in q_grid) {
      for (P in P_grid) {
        for (D in D_grid) {
          for (Q in Q_grid) {
            order <- c(p, d, q)
            seasonal <- c(P, D, Q, s)
            model <- tryCatch(
              {
                fit <- arima(west_train, order = order, seasonal = seasonal)
                AIC_val <- AIC(fit)
                # Append the results to the data frame
                west_results <- rbind(west_results, data.frame(order = paste(order, collapse = ", "),
                                                               seasonal = paste(seasonal, collapse = ", "),
                                                               AIC = AIC_val))
              },
              error = function(e) NULL
            )
          }
        }
      }
    }
  }
}

Code

save(east_results, midwest_results, south_results, west_results, file = "results/sarima_grid.rda")

Code

load("results/sarima_grid.rda")

Best East SARIMA

Code

# Find the best model (minimum AIC)
best_east_model <- east_results[which.min(east_results$AIC), ]
print(best_east_model)

      order    seasonal      AIC
213 5, 0, 4 0, 1, 1, 12 11865.53

Best Midwest SARIMA

Code

# Midwest Results
# Find the best model (minimum AIC)
best_midwest_model <- midwest_results[which.min(midwest_results$AIC), ]
print(best_midwest_model)

      order    seasonal      AIC
385 5, 1, 5 0, 0, 1, 12 12819.82

Best South SARIMA

Code

# South Results
# Find the best model (minimum AIC)
best_south_model <- south_results[which.min(south_results$AIC), ]
print(best_south_model)

      order    seasonal      AIC
402 5, 1, 5 0, 0, 0, 12 13698.93

Best West SARIMA

Code

# Find the best model (minimum AIC)
best_west_model <- west_results[which.min(west_results$AIC), ]
print(best_west_model)

      order    seasonal      AIC
385 5, 0, 5 1, 1, 0, 12 12141.63

7.4 Rebuilding SARIMA

Code

# Rebuilding SARIMA Models -----------------------------------------------------
east_sarima <- arima(east_train_ts, order = c(5, 0, 4), 
                     seasonal = list(order = c(0, 1, 1), period = 12))
midwest_sarima <- arima(midwest_train_ts, order = c(5, 1, 5),
                        seasonal = list(order = c(0, 0, 1), period = 12))
south_sarima <- arima(south_train_ts, order = c(5, 1, 5),
                           seasonal = list(order = c(0, 0, 0), period = 12))
west_sarima <- arima(west_train_ts, order = c(5, 0, 5),
                     seasonal = list(order = c(1, 1, 0), period = 12))


## Forecasting with the NEW SARIMA model ---------------------------------------
east_values <- forecast(east_sarima, h = length(east_test_ts))
midwest_values <- forecast(midwest_sarima, h = length(midwest_test_ts))
south_values <- forecast(south_sarima, h = length(south_test_ts))
west_values <- forecast(west_sarima, h = length(west_test_ts))

Code

## Plotting & Accuracy

## East Region
east_tuned_sarima_acc <- forecast::accuracy(east_values)

east_tuned_sarima_acc %>% 
  kbl(caption = "<center>Tuned SARIMA East Accuracy<center>") %>% 
  kable_styling()

Tuned SARIMA East Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.1031443	87.81972	53.94287	NaN	Inf	0.8273471	-0.0998164

Code

east_tuned_sarima_plot <- forecast::autoplot(east_values) +
  geom_line(data = east_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "East Region Tuned SARIMA",
       y = "Daily Deaths")

east_tuned_sarima_plot

Code

## Midwest Region
midwest_tuned_sarima_acc <- forecast::accuracy(midwest_values)

midwest_tuned_sarima_acc %>% 
  kbl(caption = "<center>Tuned SARIMA Midwest Accuracy<center>") %>% 
  kable_styling()

Tuned SARIMA Midwest Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	-0.1147262	143.4927	91.79687	NaN	Inf	0.4106378	-0.012698

Code

midwest_tuned_sarima_plot <- forecast::autoplot(midwest_values) +
  geom_line(data = midwest_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "Midwest Region Tuned SARIMA",
       y = "Daily Deaths")

midwest_tuned_sarima_plot

Code

## South Region
south_tuned_sarima_acc <- forecast::accuracy(south_values)

south_tuned_sarima_acc %>% 
  kbl(caption = "<center>Tuned SARIMA South Accuracy<center>") %>% 
  kable_styling()

Tuned SARIMA South Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	-4.199063	200.7617	120.224	NaN	Inf	0.3743912	-0.0666043

Code

south_tuned_sarima_plot <- forecast::autoplot(south_values) +
  geom_line(data = south_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "South Region Tuned SARIMA",
       y = "Daily Deaths")

south_tuned_sarima_plot

Code

## West Region
west_tuned_sarima_acc <- forecast::accuracy(west_values)

west_tuned_sarima_acc %>% 
  kbl(caption = "<center>Tuned SARIMA West Accuracy<center>") %>% 
  kable_styling()

Tuned SARIMA West Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	-0.0032713	108.9033	71.27948	NaN	Inf	0.4530221	0.0086743

Code

west_tuned_sarima_plot <- forecast::autoplot(west_values) +
  geom_line(data = west_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "West Region Tuned SARIMA",
       y = "Daily Deaths")

west_tuned_sarima_plot

Code

# saving results
save(east_tuned_sarima_acc, east_initial_sarima_acc, east_initial_sarima_plot, east_tuned_sarima_plot,
     midwest_tuned_sarima_acc, midwest_initial_sarima_acc, midwest_initial_sarima_plot, midwest_tuned_sarima_plot,
     south_tuned_sarima_acc, south_initial_sarima_acc, south_initial_sarima_plot, south_tuned_sarima_plot,
     west_tuned_sarima_acc, west_initial_sarima_acc, west_initial_sarima_plot, west_tuned_sarima_plot,
     file = "results/uni_sarima.rda"
     )

8. Auto ARIMA Model

8.1 Building AutoARIMA

Code

# Building AutoARIMA Model -----------------------------------------------------
east_arima <- auto.arima(east_train_ts, seasonal = TRUE)
midwest_arima <- auto.arima(midwest_train_ts, seasonal = TRUE)
south_arima <- auto.arima(south_train_ts, seasonal = TRUE)
west_arima <- auto.arima(west_train_ts, seasonal = TRUE)

# Summary of the ARIMA model
summary(east_arima)

Series: east_train_ts 
ARIMA(5,1,2) 

Coefficients:
         ar1      ar2      ar3      ar4      ar5      ma1     ma2
      0.5142  -0.7304  -0.1394  -0.3479  -0.2207  -1.0295  0.7353
s.e.  0.0480   0.0344   0.0426   0.0329   0.0418   0.0410  0.0265

sigma^2 = 6861:  log likelihood = -5956.39
AIC=11928.79   AICc=11928.93   BIC=11968.22

Training set error measures:
                    ME     RMSE     MAE MPE MAPE      MASE       ACF1
Training set 0.1506172 82.50546 47.1043 NaN  Inf 0.7224607 0.01707355

Code

summary(midwest_arima)

Series: midwest_train_ts 
ARIMA(5,0,3) with zero mean 

Coefficients:
         ar1      ar2      ar3      ar4      ar5      ma1     ma2      ma3
      0.1952  -0.6269  -0.3389  -0.2999  -0.4563  -1.2016  0.9332  -0.3022
s.e.  0.0607   0.0312   0.0356   0.0276   0.0442   0.0702  0.0718   0.0382

sigma^2 = 21712:  log likelihood = -6544.65
AIC=13107.31   AICc=13107.49   BIC=13151.67

Training set error measures:
                    ME     RMSE      MAE MPE MAPE      MASE       ACF1
Training set 0.3936445 146.7706 85.18588 NaN  Inf 0.3810646 -0.0336952

Code

summary(south_arima)

Series: south_train_ts 
ARIMA(5,0,4) with zero mean 

Coefficients:
         ar1      ar2      ar3      ar4      ar5      ma1     ma2      ma3
      0.1881  -0.9049  -0.0800  -0.4542  -0.5567  -0.9401  1.0118  -0.7340
s.e.  0.0349   0.0292   0.0464   0.0251   0.0303   0.0357  0.0434   0.0525
         ma4
      0.4897
s.e.  0.0402

sigma^2 = 40663:  log likelihood = -6864.82
AIC=13749.64   AICc=13749.86   BIC=13798.93

Training set error measures:
                    ME     RMSE      MAE MPE MAPE      MASE        ACF1
Training set 0.1422268 200.7611 119.8731 NaN  Inf 0.3732984 -0.06523796

Code

summary(west_arima)

Series: west_train_ts 
ARIMA(4,0,3) with zero mean 

Coefficients:
          ar1     ar2      ar3      ar4      ma1      ma2     ma3
      -0.0876  0.1258  -0.5602  -0.4779  -0.5985  -0.3463  0.5343
s.e.   0.0505  0.0450   0.0330   0.0309   0.0509   0.0784  0.0421

sigma^2 = 11991:  log likelihood = -6241.45
AIC=12498.9   AICc=12499.04   BIC=12538.33

Training set error measures:
                    ME     RMSE      MAE MPE MAPE      MASE        ACF1
Training set 0.1545496 109.1261 72.50105 NaN  Inf 0.4607859 -0.08454213

Code

## Forecasting Test Data -------------------------------------------------------
east_values <- forecast(east_arima, h = length(east_test_ts))
midwest_values <- forecast(midwest_arima, h = length(midwest_test_ts))
south_values <- forecast(south_arima, h = length(south_test_ts))
west_values <- forecast(west_arima, h = length(west_test_ts))

8.2 Accuracy

Code

## Plotting & Accuracy

## East Region
east_auto_acc <- forecast::accuracy(east_values)

east_auto_acc  %>% 
  kbl(caption = "<center>Auto ARIMA East Accuracy<center>") %>% 
  kable_styling()

Auto ARIMA East Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.1506172	82.50546	47.1043	NaN	Inf	0.7224607	0.0170736

Code

east_auto_plot <- forecast::autoplot(east_values) +
  geom_line(data = east_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "East Region Auto ARIMA",
       y = "Daily Deaths")

east_auto_plot 

Code

## Midwest Region
midwest_auto_acc <- forecast::accuracy(midwest_values)

midwest_auto_acc %>% 
  kbl(caption = "<center>Auto ARIMA Midwest Accuracy<center>") %>% 
  kable_styling()

Auto ARIMA Midwest Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.3936445	146.7706	85.18588	NaN	Inf	0.3810646	-0.0336952

Code

midwest_auto_plot <- forecast::autoplot(midwest_values) +
  geom_line(data = midwest_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "Midwest Region Auto ARIMA",
       y = "Daily Deaths")

midwest_auto_plot 

Code

## South Region
south_auto_acc <- forecast::accuracy(south_values)

south_auto_acc %>% 
  kbl(caption = "<center>Auto ARIMA South Accuracy<center>") %>% 
  kable_styling()

Auto ARIMA South Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.1422268	200.7611	119.8731	NaN	Inf	0.3732984	-0.065238

Code

south_auto_plot <- forecast::autoplot(south_values) +
  geom_line(data = south_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "South Region Auto ARIMA",
       y = "Daily Deaths")

south_auto_plot

Code

## West Region
west_auto_acc <- forecast::accuracy(west_values)

west_auto_acc %>% 
  kbl(caption = "<center>Auto ARIMA West Accuracy<center>") %>% 
  kable_styling()

Auto ARIMA West Accuracy

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	0.1545496	109.1261	72.50105	NaN	Inf	0.4607859	-0.0845421

Code

west_auto_plot <- forecast::autoplot(west_values) +
  geom_line(data = west_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "West Region Auto ARIMA",
       y = "Daily Deaths")

west_auto_plot

Code

# saving results
save(east_auto_plot,
     midwest_auto_plot,
     south_auto_plot,
     west_auto_plot,
     east_auto_acc,
     midwest_auto_acc,
     south_auto_acc,
     west_auto_acc,
     file = "results/uni_autoarima.rda")

9. Compare Arima Models

Code

load("results/uni_arima.rda")
load("results/uni_sarima.rda")
load("results/uni_autoarima.rda")

East Results

Code

east1 <- east_initial_arima_acc %>% 
  as.data.frame()

east2 <- east_tuned_arima_acc %>% 
  as.data.frame()

east3 <- east_initial_sarima_acc%>% 
  as.data.frame() 

east4 <- east_tuned_sarima_acc%>% 
  as.data.frame()

east5 <- east_auto_acc%>% 
  as.data.frame()

tribble(
  ~model, ~mae, ~mase,
  #--|--|----
  "initial arima", east1$MAE, east1$MASE,
  "tuned arima", east2$MAE, east2$MASE,
  "initial sarima", east3$MAE, east3$MASE,
  "tuned sarima", east4$MAE, east4$MASE,
  "auto arima", east5$MAE, east5$MASE,
) %>% 
  kbl() %>% 
  kable_styling()

model	mae	mase
initial arima	64.66739	0.9918341
tuned arima	47.18823	0.7237480
initial sarima	61.27670	0.9398295
tuned sarima	53.94287	0.8273471
auto arima	47.10430	0.7224607

Midwest Results

Code

midwest1 <- midwest_initial_arima_acc %>% 
  as.data.frame()

midwest2 <- midwest_tuned_arima_acc %>% 
  as.data.frame()

midwest3 <- midwest_initial_sarima_acc%>% 
  as.data.frame() 

midwest4 <- midwest_tuned_sarima_acc%>% 
  as.data.frame()

midwest5 <- midwest_auto_acc%>% 
  as.data.frame()

tribble(
  ~model, ~mae, ~mase,
  #--|--|----
  "initial arima", midwest1$MAE, midwest1$MASE,
  "tuned arima", midwest2$MAE, midwest2$MASE,
  "initial sarima", midwest3$MAE, midwest3$MASE,
  "tuned sarima", midwest4$MAE, midwest4$MASE,
  "auto arima", midwest5$MAE, midwest5$MASE,
) %>% 
  kbl() %>% 
  kable_styling()

model	mae	mase
initial arima	138.93519	0.6215031
tuned arima	95.00769	0.4250009
initial sarima	138.93519	0.6215031
tuned sarima	91.79687	0.4106378
auto arima	85.18588	0.3810646

South Results

Code

south1 <- south_initial_arima_acc %>% 
  as.data.frame()

south2 <- south_tuned_arima_acc %>% 
  as.data.frame()

south3 <- south_initial_sarima_acc%>% 
  as.data.frame() 

south4 <- south_tuned_sarima_acc%>% 
  as.data.frame()

south5 <- south_auto_acc%>% 
  as.data.frame()

tribble(
  ~model, ~mae, ~mase,
  #--|--|----
  "initial arima", south1$MAE, south1$MASE,
  "tuned arima", south2$MAE, south2$MASE,
  "initial sarima", south3$MAE, south3$MASE,
  "tuned sarima", south4$MAE, south4$MASE,
  "auto arima", south5$MAE, south5$MASE,
) %>% 
  kbl() %>% 
  kable_styling()

model	mae	mase
initial arima	203.7282	0.6344327
tuned arima	120.2240	0.3743912
initial sarima	203.7282	0.6344327
tuned sarima	120.2240	0.3743912
auto arima	119.8731	0.3732984

West Results

Code

west1 <- west_initial_arima_acc %>% 
  as.data.frame()

west2 <- west_tuned_arima_acc %>% 
  as.data.frame()

west3 <- west_initial_sarima_acc%>% 
  as.data.frame() 

west4 <- west_tuned_sarima_acc%>% 
  as.data.frame()

west5 <- west_auto_acc%>% 
  as.data.frame()

tribble(
  ~model, ~mae, ~mase,
  #--|--|----
  "initial arima", west1$MAE, west1$MASE,
  "tuned arima", west2$MAE, west2$MASE,
  "initial sarima", west3$MAE, west3$MASE,
  "tuned sarima", west4$MAE, west4$MASE,
  "auto arima", west5$MAE, west5$MASE,
) %>% 
  kbl() %>% 
  kable_styling()

model	mae	mase
initial arima	101.08083	0.6424269
tuned arima	63.93076	0.4063168
initial sarima	101.08083	0.6424269
tuned sarima	71.27948	0.4530221
auto arima	72.50105	0.4607859